Shell-Model Monte Carlo Simulations of BCS-BEC Crossover in Few-Fermion Systems 
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We study a trapped system of fermions with a zero-range two-body interaction using the shell- 
model Monte Carlo method, providing ab initio results for the low particle number limit where 
mean-field theory is not applicable. We present results for the iV-body energies as function of 
interaction strength, particle number, and temperature. The subtle question of renormalization in 
a finite model space is addressed and the convergence of our method and its applicability across 
the BCS-BEC crossover is discussed. Our findings indicate that very good quantitative results can 
be obtained on the BCS side, whereas at unitarity and in the BEC regime the convergence is less 
clear. Comparison to N — 2 analytics at zero and finite temperature, and to other calculations in 
the literature for N > 2 show very good agreement. 

PACS numbers: 03.75.Ss, 05.30.Fk, 21.60.Ka 



I. INTRODUCTION 

The physics of ultracold atomic gases has been inten- 
sively pursued experimentally and theoretically in the 
last decade. Recently there has been great interest in 
strongly-interacting Fermi gases where Feshbach reso- 
nances allow the tuning of the two-body interaction, and 
studies of the transition from a dilute gas of fermionic 
atoms to a Bose condensate of molecules are now possi- 
ble in the laboratory 0, H i, H H 0. While studies of 
degenerate Fermi gases have mostly dealt with large atom 
numbers and wide traps, efforts have begun to trap only 
a few atoms (1-100) in tighter traps [7|. Also, with the 
implementation of three-dimensional optical lattices, a 
low-tunneling regime can be reached with essentially iso- 
lated harmonic oscillators containing only a few fermions 
at each site Q. This means that one can now explore 
few-body fermionic effects in trapped systems with scat- 
tering lengths that are comparable to the inter-particle 
distance and the trap width. 

In this paper we report on a theoretical study of har- 
monically trapped fermions using the shell-model Monte 
Carlo (SMMC) approach. This method has been ex- 
tensively used in nuclear physics to determine nuclear 
properties at finite temperature in larger model spaces 
than can be handled by normal nuclear shell-model diag- 
onalization @, [l(J. In the SMMC, the many-body prob- 
lem is described by a canonical ensemble at tempera- 
ture IcbT — f3 and the Hubbard-Stratonovich trans- 
formation is used to linearize the imaginary-time many- 



body propagator e~@ H . Observables are then expressed 
as path integrals of one-body propagators in fluctuating 
auxiliary fields. The method is in principle exact and 
subject only to statistical uncertainties. For equal mix- 
tures of two hyperfine states at low density, the inter- 
action can be modeled with an s-wave zero-range poten- 
tial. Importantly, this interaction is free of sign problems 
[ll( that are otherwise known to plague quantum Monte 
Carlo simulations with fermions. We present here the 
first application of this many-body method to ultracold 
gas physics. 
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Previous works have considered few-fermion systems 
using advanced many-body methods. The Green's func- 
tion Monte Carlo methods were applied to homogeneous 
[13], as well as trapped systems [la ]. No-core [3], and 
traditional shell- models [HI, using effective interactions 
have also recently been applied to these systems, particu- 
larly for very low particle numbers where exact results are 
available [la. fT?|. Finite-temperature, non-perturbative 
lattice methods have also been applied to the homoge- 
neous case fl9j . These works mostly focus on the 
unitary |a| — > oo limit and the crossover regime around it. 
Most previous Monte Carlo approaches have used fixed 
nodes in the many-body wave function in order to alle- 
viate the sign problem, making the methods variational. 
As we will now demonstrate the present method has no 
sign problem [ll[ and can be used on the BCS side, in 
the crossover region, and also into the BEC regime. 
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II. MODEL AND RENORMALIZATION 
SCHEME 

The model Hamiltonian used is 

H = E f| + E l m ^ r i + E VoSft r 3 ), (1) 
i to] 

where we sum over all particles i and [zj] denotes a 
sum over fermion pairs with opposite internal (hyper- 
fine) states. The trap frequency is lo, and Vo denotes the 
interaction strength. The SMMC method was originally 
set up to handle nucleons where the Hamiltonian above 
appears in extensively studied pairing problems in nuclei 
[20( . The two-component Fermi gas can now be mapped 
onto a single spin 1/2 nucleon species [TTI. [2l| . 

Dimensional arguments reveal that the matrix ele- 
ments needed in a shell-model approach scale as 1/b 3 , 
where b = ^jTi/mui is the oscillator length. It is there- 
fore natural to redefine the interaction strength in terms 
of Vo — —ghuib 3 , where g is a dimensionless strength 
measure. In order to relate to physical quantities the 
strength of the zero-range interaction must be regularized 
[22| . The shell-model works in finite model spaces and 
this naturally introduces a cut-off in energy E c = a 2 hu>, 
where a 2 = N max + 3/2. Regularization in finite model 
spaces with discrete energies is notoriously difficult and 
several prescriptions have been adopted in the literature. 
Here we will use the simplest strategy and renormalize 
the coupling g through low-energy scattering parameters 
defined in the continuum [23[ . This defines a relation be- 
tween g 7 the s-wave scattering length a, and E c , which 
is Ana/b = g/(ag/2ir 2 — 1). This yields the effective in- 
teraction strength for a model space with N max + 1 ma- 
jor harmonic oscillator shells and gives a prescription for 
varying the interaction strength with model space size in 
order to keep a fixed. 

The relevant parameter regime is expected to be where 
the natural energy scale, given by the level spacing Tito, is 
comparable to typical two-body matrix elements between 
the trap states at T = (as a quantitative measure we 
use (lsls|V|lsls)). This turns out to be around g ~ 10 
in our setup at N max = 3, which corresponds to a = 
lib in the continuum regularization scheme. The regions 
of interest in terms of interaction strength and hu> for 
atomic gases are discussed in [24|, [25| . 

It is appropriate to discuss alternative renormalization 
schemes that were recently proposed [3, EH which take 
the trap level spacing into account in a direct way in- 
spired by effective field theory techniques. The idea is to 
relate the bare coupling Vq to a by calculating the two- 
body ground state energy in the finite model space of 
harmonic oscillator states with energy less than E c and 
then determining a through the exactly solvable pseu- 
dopotential model of Busch et al. (2(|, i.e. solving the 
equation T(3/4:-E/2u))/T(l/4-E/2w) = b/V^a, with T 
the Gamma function. There has been some experimental 
support of the Busch et al. results for fermion pairs in 
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FIG. 1: (Color online) Running of the coupling, — ag/4ir, 
in continuum and pseudopotential regularization schemes as 
function of the inverse cut-off (v^a)- 1 . The upper panel 
shows the a/b = —1.0 (BCS) and a/b = 11 (BEC) cases used 
in the calculations. The dashed line is the continuum scheme 
with the couplings g = 10 (BEC) and g = 5.35 (BCS) at 
N max = 3, whereas the solid line is the pseudopotential line 
that runs through the same couplings at N max = 3. The 
dotted line is the pseudopotential coupling for the same a/b 
ratio as in the continuum case. The dot-dashed line is uni- 
tarity \a\ = oo. The lower panel demonstrates the deviations 
in the deep BEC regime a/b = 1.0. Notice the difference in 
scales on the vertical axis between the upper and lower panels. 



three-dimensional (3D) optical lattices 8]. However, the 
pseudopotential is only expected to be correct when the 
van der Waals scale is much smaller than the trap, and 
since we have tight traps with only few fermions in mind 
this is not necessarily the case. We therefore expect both 
the continuum and the pseudopotential approaches to be 
approximations for a tight trap [27] . In Fig. [1] we have 
plotted the running of the coupling in both schemes simi- 
larly to Fig. 1 of The upper part shows the two cases 
used in our calculations, whereas the lower part is in the 
deep BEC regime (a/b = 1.0). The two approaches yield 
a/b = 11 (continuum) and a/b = 7.8 (pseudopotential) 
for g — 10 with N max = 3, which is not too severe con- 
sidering the divergence of a at unitarity. In the case of 
g = 5.35 at N max = 3 we have a/b = —1.0 (continuum) 
and a/b = —0.9 (pseudopotential). In the exact model 
of [201 this translates to a difference in ground state en- 
ergy of only 0.03K.W for both cases. The running of the 
coupling is the same to two significant digits for both 
a/b = 11 and a/b = —1.0 results for N max — 2, 3, and 4 
used here. If we turn the argument around and instead 
fix the ratio a/b, the pseudopotential couplings are given 
as the dashed blue curves in Fig. [T] For the a/b = —1.0 
case, this gives g = 5.67 as compared to our g = 5.35, 
and g = 9.79 compared to g — 10 for a/b — 11. These 
modifications would only bring us slightly downward on 
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TABLE I: Energies (in units of hui) calculated with the 
SMMC method for a trapped fermion gas with scattering 
lengths a/b = 11 (BEC) and a/b = -1.0 (BCS) calculated 
at temperature ksT = l/5hu for different particle numbers 
TV. The statistical uncertainty is given in parenthesis. HOSD 
denotes the non- interacting energies at T = 0. 



N 


HOSD 


BEC BCS 


N 


HOSD 


BEC 


BCS 


2 


3 


1.72(3) 2.49(3) 


12 


32 


20.7(2) 


27.23(2) 


3 


5.5 
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FIG. 2: (Color online) Scaled energies E mt /N 4 ' 3 , where 
E int = E - E g=0 (upper panel), and A(7V) = E(N + 1) - 
2E(N) + E(N — 1) (lower panel) as a function of particle 
number N for a/b — 11 (solid) and a/b — —1.0 (dashed). 



the BCS side and slightly upward close to unitarity. Only 
on the deeper BEC side of the crossover (5 > a/b > in 
the continuum scheme), where the pairs are essentially 
molecules, do we find larger deviations of the coupling in 
the two schemes. The lower panel in Fig. [T] shows the 
a/b = 1.0 case. Here one can clearly see that the run- 
ning of the schemes is very different expect for very large 
model spaces. We do not expect the SMMC approach to 
work in the molecular regime. 



III. RESULTS 
A. Zero-temperature Limit 

In TableUwe present the SMMC energies for two values 
of the scattering length a, along with the non-interacting 
energies, for particle numbers N = 2 — 21. All results 
in the table were calculated with a cut-off N max = 3, 
corresponding to a four major shell model space. The 
results at low T are well converged in this model space 
for a/b = —1.0 as we demonstrate below. For a/b = 11 
the energies suffer larger uncertainties as discussed be- 
low. The results in the table were all calculated at in- 
verse temperature [3 = 5/fku, above which we find only 
statistical changes in energy (see below). Thus the re- 
sults in Table U are the SMMC estimates of the T = 
energy. To explore the results, we plot (E — E g= o)/N 4 / 3 
vs. N in the upper part of Fig. [51 which is the in- 
teraction energy divided by the Thomas-Fermi scaling. 
We can clearly see an odd-even staggering that becomes 
more pronounced on the BEC side. The lower panel in 
Fig. [5] shows this more clearly through the fundamental 
gap A(N) = E(N + 1)-2E(N)+E(N-1). Notice here 
that the closed shells are prominent and even more so 
on the BCS side. We thus predict that energy measure- 
ments alone can probe the crossover through odd-even 



staggering at and away from closed shells. 

As discussed above, the BEC results in Table H] were 
obtained with g = 10 at N max = 3 or a = lib which 
is on the BEC side of the crossover but still quite close 
to unitarity which is located at g = 9.31 for N max — 3 
(for both regularization schemes discussed above). The 
energies at unitarity are therefore slightly above our BEC 
results. The BCS results were obtained with g = 5.346 
or a = —1.06. 

As discussed, the N = 2 problem can be solved exactly 
for all a in a pseudopotential approach (2(| . If we insert 
a = lib we find E — 1.92huj for the T — ground state, 
whereas for a = —1.06 we find E — 2A9frui. Our values of 
1.72(3) and 2.49(3) at T = l/hhu/k B thus indicates that 
we are close to the T = limit. We notice that the energy 
is underestimated on the BEC side, indicating that our 
method becomes worse in that regime. 

In general we find excellent agreement with previous 
calculations. For N = 3 the exact energy at unitarity 
is E = A.27hw [l6j]. The results in Table U are slightly 
below (above) this value on the BEC (BCS) side, as we 
would expect. For higher N, we find good agreement 
with the results of [28] . A comparison can also be made 
with the results at unitarity of [13j which cover the par- 
ticle numbers presented here, and for which our a = lib 
energies are slightly smaller but sandwiched between the 
values in Table II of [HI]. In order to elaborate on the 
model space size effects we plot in Fig. [3] the energy as 
function of N max for N — 2, 3, and 4 at the converged 
temperature (3 = huj/5 (see below) for a/b = 11 (upper 
panel) and a/b = —1.0 (lower panel). The convergence is 
clearly much better on the BCS than on the BEC side. 
In particular, we see a decreasing trend of the energy 
with model space size for a/b = 11, although we remark 
that at least for N = 2 and 3, the results are consistent 



4 



3 4 

6q 3 
























. 5 - a/b=-1.0 



N„ 



FIG. 3: (Color online) Energy in units of Twj as function 
of model space size N max as denned in the text for parti- 
cle numbers N = 2, 3, and 4 with a/b = 11 (upper panel) and 
a/b — —1.0 (lower panel). 
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FIG. 4: (Color online) Energies calculated with the SMMC 
method in units of hui as a function of inverse temperature 
/3 = (fcsT) -1 for N = 2, 6, and 10 with strength parameter 
a/b = 11 (upper panel) and a/b = —1.0 (lower panel) for 
model space sizes N max — 2 (black), 3 (red), and 4 (blue). 



within the statistical uncertainties. We thus see that the 
SMMC is robust in calculating the ground state energies 
of small Fermi systems on the BCS side. At unitarity 
and into the BEC regime our results are not as robust, 
and the deep BEC regime cannot be accessed. 

B. Model space size and convergence 

The SMMC method works at finite temperature and 
T = results are obtained by increasing (3 = (kgT) -1 
and finding the point of convergence. As we work in fi- 
nite model spaces and must use regularized interactions, 
it is imperative to have this convergence under control. 
However, as the SMMC is also a shell-model method, 
there is a fundamental computational restriction on the 
model space size that can be used. The SMMC requires 
diagonalization of one-body terms, which grow quadrat- 
ically in the number of model space states. While being 
more benign in its scaling, a compromise must still be 
found that allows many calculations at different (3 val- 
ues with good Monte Carlo statistics. We caution that 
larger model spaces also lead to more numerical noise in 
the calculation, and that this will grow as the coupling 
gets stronger. 

Figure 2] shows the SMMC energies as a function of 
(3 for selected even particle numbers N = 2, 6, and 10 
for a = lib (upper panel) and a = —1.06 (lower panel). 
The horizontal bars indicate the error estimates from the 
Monte Carlo integrations. The first thing we notice is 
that the energies are virtually converged at (3 > 5/hoj 
for all N and all model spaces in the a/b = —1.0 case. 
The fluctuations above this point are largely due to nu- 
merical noise at large (3 that makes sampling difficult. 
This justifies our choice of (3 = 5/huj in Table U We 
see the convergence getting worse with N, which is ex- 



pected as there will be model space saturation. However, 
within the statistical uncertainties, we see good agree- 
ment for the different model spaces for both N = 2 and 
6. For N ~ 10 there is a slight decrease of the energy 
with model space size and our convergence is not quite as 
good. This is clear since for larger N there is less space 
to excite particles in the given model space. For a/b = 11 
the convergence is noticeably worse with larger fluctua- 
tions. Here the pairing is strong and a mixing of scales 
in the numerics makes the Monte Carlo sampling diffi- 
cult. This problem becomes worse for larger g. Relating 
back to Table U the N > 10 results are therefore proba- 
bly slightly overestimated. However, the good agreement 
with the results of [13[ indicates that the deviations are 
under control and on the scale of that seen for N = 10. 

For odd N the projection onto exact number states 
introduces a sign problem [l(| which grows with the cou- 
pling. Therefore the odd N results in Table U tend to 
have larger uncertainties. This can also be seen in Fig. [3] 
Notice that the large differences at low (3 in Fig. Q] are 
due to the finite model space. At high T, the particles 
will equilibrate in the available states and the energy be- 
comes E ~ NE, where E — J2i e i/di with e, the zth 
single-particle energy and di the corresponding degener- 
acy. This will in turn make the energy in the high-T 
limit grow with N max toward its thermodynamic value 
E(T) = SNksT for the non-interacting harmonic oscil- 
lator. 

Throughout the calculation we have used the contin- 
uum regularization method. As noted, only on the deep 
BEC side will there be large differences to the pseu- 
dopotential regularization approach of [3, EB] ■ However, 
these renormalization schemes do not depend explicitly 
on JV. In our calculations we see the convergence be- 
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FIG. 5: (Color online) Energy for N = 2 as function of tem- 
perature, T, in units of hu for a/6 = 11 (left) and a/6 = —1.0 
(right) with model spaces N max = 2 (black), N max = 3 (red), 
and N m ax = 4 (blue). The black dashed line is the Busch et 
al. [U result. 



coming worse for high N which might suggest an N de- 
pendence; Vq = Vq (a, b, N, N max ) . This can be improved 
upon by using effective interactions in finite model spaces 
as done often in nuclear physics where good effective two- 
body interactions do depend on N. We will explore such 
options in the future. 

C. Finite temperature N = 2 sesults 

To explore the temperature dependence of our results, 
we present in Fig. [5] the N — 2 energies with a/6 = 11 
(left) and a/6 = —1.0 (right) for different model spaces. 
The dashed line in the figure is obtained from the Busch 
et al. [2(| results. Here we see that the SMMC obtains 
very good agreement as function of T with the exact re- 
sult, except for the very low and high T regions. At 
low T we see deviations and fluctuations for a/b = 11. 



The latter are due to numerical problems at strong cou- 
pling, whereas we believe the deviation from the exact 
result of the average is connected to the regularization as 
discussed above. For intermediate T we find very good 
agreement with the exact results, and as in Fig. [4] we 
also see that we have obtained convergence in our larger 
model spaces below T ~ 0.8huj/kB ■ The latter is very 
important for finite-temperature quantities such as the 
specific heat, where model space sizes can produce so- 
called Schottky peaks [Io| . For N > 2 we have indica- 
tions of a pairing phase transition in the pair correlations 
and the specific heat. These results will be presented 
shortly. 

IV. CONCLUSION AND OUTLOOK 

As we have shown, the SMMC offers a good quantita- 
tive description of the behavior of small Fermi systems in 
an interesting interaction regime. We studied the case of 
isotropic traps and calculated the many-body energies. 
The excellent convergence properties of the method on 
the BCS side of the BCS-BEC crossover holds promise 
for application to specific situations where, e.g., defor- 
mation properties and formation of higher angular mo- 
mentum pairs may become relevant. Deformation was 
studied in the nuclear case using the SMMC to predict 
shape transitions occurring as a function of temperature 
in the competition between pairing and quadrupole inter- 
actions [10| . This is relevant also for ultracold gases with 
the recent realization of condensates with intrinsic long- 
range interactions between the atoms [29(. Recently the 
SMMC was also used to study the parity and spinprop- 
erties of the density of states in nuclear systems [30, [31| . 
Transforming this to the atomic system could help us 
understand the low-energy excitation spectrum and the 
response of the gas to external perturbations. 
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